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ABSTRACT 

Galactic cosmic rays are believed to be accelerated at supernova remnant shocks. 
Though very popular and robust, this conjecture still needs a conclusive proof. The 
strongest support to this idea is probably the fact that supernova remnants are ob- 
served in gamma-rays, which are indeed expected as the result of the hadronic in- 
teractions between the cosmic rays accelerated at the shock and the ambient gas. 
However, also leptonic processes can, in most cases, explain the observed gamma-ray 
emission. This implies that the detections in gamma rays do not necessarily mean 
that supernova remnants accelerate cosmic ray protons. To overcome this degeneracy, 
the multi-wavelength emission (from radio to gamma rays) from individual super- 
nova remnants has been studied and in a few cases it has been possible to ascribe 
the gamma-ray emission to one of the two processes (hadronic or leptonic) . Here we 
adopt a different approach and, instead of a case-by-case study we aim for a popu- 
lation study and we compute the number of supernova remnants which are expected 
to be seen in TeV gamma rays above a given flux under the assumption that these 
objects indeed are the sources of cosmic rays. The predictions found here match well 
with current observational results, thus providing a novel consistency check for the 
supernova remnant paradigm for the origin of galactic cosmic rays. Moreover, hints 
are presented for the fact that particle spectra significantly steeper than E~ 2 are 
produced at supernova remnants. Finally, we expect that several of the supernova 
remnants detected by H.E.S.S. in the survey of the galactic plane should exhibit a 
gamma-ray emission dominated by hadronic processes (i.e. neutral pion decay). The 
fraction of the detected remnants for which the leptonic emission dominates over the 
hadronic one depends on the assumed values of the physical parameters (especially 
the magnetic field strength at the shock) and can be as high as roughly a half. 
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1 INTRODUCTION 



Though cosm ic Rays (CRs) ha ve been discovered about a 
century ago (|De Angelid l2012h . the question about their 



origin is still a matter of d iscussion (see e.g. Drurv et al 



l200ll; [Ah aronian et al.ll2012l . for reviews). iBaade fc Zwickvl 



( 1934) have been the first to propose that supernovae are the 
sources of CRs, and this still remains the most popular sce- 
nario to explain the origin of galactic CRs. The particle en- 
ergy that marks the transition between galactic and extra- 
galactic CRs is believed to be located between the knee and 
the ankle which are observed in the CR spectrum at par- 
ticle energies of « 4 PeV and a few 10 18 eV, respectively. 
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Thus, the acceleration mechanism connected to supernovae 
must be able to accelerate particles up to the PeV energy 
range and above. The present formulation of this idea is of- 
ten referred to as the SuperNova Remnant (SNR) paradigm 
for the origin of CRs, because galactic CRs are believed to 
be accelerated vi a first order Fermi mechanism o perating at 
SNR shocks (e.g. lHillasj|2005l : iHelder et al]|2012r i. 

The success of this paradigm relies on several facts. 
First of all, SNRs can provide the power needed to sustain 
the CR flux at the observed level, if some fraction (about 
10. ..30%) of their kinetic energy is somewhat converted into 
relativistic particles (e.g. iPtuskin et alj|20icl . and references 
therein). Second, diffusive shoc k acceleration can operate at 
SNR shocks (for a review see lDrurvlll983h . and this pro- 
vides a viable mechanism to accelerate CRs. Diffusive shock 
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acceleration is expected to accelerate particles with power 
law spectra whose slope, once the effects of CR propagation 
in the Galaxy are taken into account, roughly matches the 
one of the CR spectrum as it is observed at the E arth (e.g. 
iPtuskin et al. I l20ld ; iBell et alJlioIH ; ICapriolillioij ) . Finally, 
during the acceleration process CRs can amplify the mag- 
netic field at shocks via various plas ma instabilities fe.g. lBelll 
1 19781 . 120041 ; iDrurv fc Downesll2012r ) . Observa tional evide nce 
for such an amplification has been found fsee lVinkll2012l . for 
a review) , and it is likely that the magnetic field strength at 
SNR shocks may grow up to a level that allows the accel- 
eration of particles up to PeV energies and more. All these 
things support the idea that SNRs indeed are the sources 
of CRs, but it has to be kept in mind that an unambiguous 
and conclusive proof of such a statement is still missing. 

The acceleration of CRs must be accompanied by the 
production of gamma-rays. This radiation is the result of 
the decay of neutral pions generated in hadroni c inter- 
actions between the CRs and the a mbient gas (ISteckerl 
1 1971); lAharonianll200^ It was shown bv lDrurv et all (ll994T l 
and iNaito fc Takaharal l|l994 l that if SNRs indeed are the 
sources of CRs, then their gamma-ray emission must be 
strong enough to be detected by Cherenkov telescopes of 
current generation. The detection of several SNRs in TeV 
gamma rays nicely fits with these earlier predictions, but 
it cannot be considered a proof of the fact that SNRs 
can accelerate CRs. This is because electrons can as well 
be accelerated at shocks, and their inverse Compton emis- 
sion can also account for the observed TeV radiation (e.g. 
Aharonian et all l2008dl ; iGabicil 120081 : iHinton fc Hofmannl 



20091 ). The ambiguity between the hadronic or leptonic ori- 



gin of the gamma-ray emission observed from SNRs is the 
main obstacle in proving (or falsifying) the SNR paradigm 
for the origin of CRs. 

Multi-wavelength observations of SNRs, from the radio 
band to the very high energy gamma-ray domain, can help 
in solving such degeneracy. This is generally done on a case 
by case basis, i.e. multi-wavelength data are collected for a 
specific SNR, and ha dronic and/or leptonic models are fit- 
ted to data (see e.g. lEllison et al.l |2012|; iBerezhko fc Volkl 



l20ld ; IZirakashvili fc Aharonianl |2010| ; iMorlino et all 120091 . 

and references therein). In some cases, it has been possible 
to ascribe quite confidently the gamma -ray emission either 
to a hadronic or a leptoni c mechanism ( Acciari et al. 20 111: 



Morl ino fc Capriolil l2012i ; lAbdo et all 



201 ll ; lEllison et al 



2010D . while for other cases the situation still remains am- 



biguous. 

In this paper, we follow a different approach and, in- 
stead of considering one specific object, we investigate the 
gamma-ray properties of SNRs as a class of objects. We start 
from the assumption that SNRs are the sources of CRs. This 
assumption, together with the knowledge of the supernova 
rate in the Galaxy, allows us to infer the typical CR acceler- 
ation efficiency per SNR. Then, by means of a Monte Carlo 
method we simulate the location and time of explosion of 
all the supernovae in the Galaxy. Finally, from the informa- 
tion on the gas density, taken from galactic surveys of CO 
and HI lines (that trace molecular and atomic Hydrogen, re- 
spectively), it is possible to estimate the hadronic gamma- 
ray emission from each simulated SNR. A leptonic contri- 
bution is added, by treating the electron-to-proton ratio 
K ep as a free parameter of the model. Following this proce- 



dure it is possible to build mock-catalogues of TeV-bright 
SNRs that can then be compared with the data coming, 
for example, from the H.E.S.S. survey of the galactic plane 
|Ah aroman et al.| [2005V Our results show that expectations 
match quite well current observations, providing an addi- 
tional and novel consistency check for the SNR paradigm 
for the origin of galactic cosmic rays. 

The paper is structured as follows: in Section [2] we de- 
scribe a procedure to estimate the gamma-ray emission from 
an individual SNR at a given stage of evolution, while in 
Section [3] a Monte Carlo approach is adopted to simulate 
the position and time of explosion of all the supernovae 
that exploded in the Galaxy. These results are then used 
in Section [3] to estimate the number of SNR detectable in 
the Galaxy at a given gamma-ray flux. In Section [4] a com- 
parison with existing data (mainly from the survey of the 
galactic plane performed by the H.E.S.S. collaboration) is 
performed. Finally, we discuss the results and conclude in 
Sections [S] and [U 



2 A MODEL FOR COSMIC RAY 

ACCELERATION AND GAMMA RAY 
PRODUCTION IN SUPERNOVA 
REMNANTS 

In this section we develop a model that couples the dynam- 
ical evolution of a SNR with the particle acceleration oper- 
ating at the shock. The aim of the model is to obtain pre- 
dictions for the gamma-ray emission from individual SNRs. 



2.1 Dynamical evolution of supernova remnants 

In order to determine the time evolution of the SNR shock 
radius and velocity, we follow th e approach outlined in 
IPtuskin fc Zirakashvilil (|200ll2005l ). where a significant con- 
tribution of CRs to the pressure behind the SNR shock has 
been assumed. 

Let us consider first the case of a thermonuclear, type la, 
supernova. The time evolution of the shock radius R„h and 
velocity u s h in the ejecta-dominated pha se are described 
by the following self-simila r expressions (Chevalier Il982l : 
IPtuskin fc Zirakashvilil 120051 ): 



Rsh = 5.3 



u 3h = 3.0 x 10 J 



f 2 



n M, 



ej,0 



1/7 



,4/7 
''kyr 



pc 



n M, 



1/7 

t 



-3/7 
kyr 



km/s 



(1) 



where £51 is the supernova explosion energy in units 
of 10 51 erg, no is the ambient gas number density in 
cm -3 , Me-,,0 is the mass ejected in the explosion in so- 
lar mass units, and t^ yT is the time after explosion ex- 
pressed in kilo-years. To describe the SNR evolution during 
the adiabatic phase it is convenient to use the expres sions 
dTruelove fc McKeelll999l ; IPtuskin fc Zirakashvilill2005h : 
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u sh = 1.7 x 1(T 
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which connect smoothly with Equations [T] at a time to ~ 
260(Me 7 -.p)/1.4) 5/6 g,~ 1/: ^n^ 1/3 yr and tend to the exact 
Sedov-Taylor solution (|Sedovlll959l ; lTavloilll95Ch for t > to. 
We follow the SNR evolution until the time t ra d ~ 3.6 x 
i n r] 4 ^ 7 yr, whi c h ma rks the transition to the radia- 



tive phase lcioffi et all l| 19881 ). 

Different expressions need to be adopted to describe the 
evolution of a core-collapse supernova, whose shock prop- 
agates in the wind blown bu bble generated by the wind o f 
the progenitor star. Following lPtuskin fc Zirakashvilil (|2005l ) 
we divide the wind blown bubble into two regions: a dense 
red-supergiant wind and a tenuous hot bubble which has 
been inflated by the wind of the massive progenitor star 
in main sequence. The red-supergiant wind is assumed to 
be spherically symmetric with velocity u w = W 6 u Wt e cm/s, 
mass loss rate M = 10~ M_sM©/yr, and density n w = 
M I 'A-Rm a u m r 2 , where m p is the proton mass and m a ~ pan p 
is the mean interstellar atom mass (here we adopt fi ~ 1.4) 
and r is the distance from the star. The radius of the 
wind is fixed to R w w 2 pc, since its exact location does 
not affect significantly the results. The radius of the hot 
bubble is Rb = 28(1,36 / M n o) 1//5 £Myr P c > where Lse is the 
main-sequence star wind power in units of 10 36 erg/s, and 
tMyr is the wind lifetime in units of mega-years. The den- 
sity inside the bubb le is n b = 0.0l(L^ 6 nl 9 t^ y 2 r ) 1/:i5 cm" 3 
I Weaver et al.lll977f ). Here we assume tMyr to be of the or- 
der of several Myr, which corresponds to th e duration of t he 
main sequence phase of very massive stars (Longair 1 201 if ). 

In such a structured interstellar medium, the evolution 
of the SNR shock during the ejecta dominated phase is de- 
scribed by (|Chevalieilll982l ; iPtuskin fc Zirakashvilill2005l ): 
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Fairly accurate expressions that describe the evolution 
of a SNR in the adiabatic phase can be obtained, in 
this case, by adopting t he thin-shell approximation (e 
lOstriker fc McKed fl995l : iBisnovatvi-Kogan fc Silichlll99, 
i.e. the assumption that the gas swept up by the SNR shock 
is concentrated in a thin layer behind the shock front. The 
following equations can be derived, where the shock speed 
and S NR age are parametrized as functions of the shock 
radius (jPtuskin fc Zirakashvilill2005l ): 



Ush(Rsh) 



7ad + 1 
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6(7ad-l)/(7ad + l) 
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Jo u sh {r) 

In the expressions above, 7 a( j is the gas adiabatic index 



and M(R s h) is the total gas mass inside the SNR shock 
(ejecta+swept-up). Equations [3] and [4] are fitted together at 
the transition between the ejecta-dominated and adiabatic 
phases. We follow the evolution of the SNR until the Mach 
number drops to ~ 3, or, if shorter, until the time at which 
the SNR shock impacts onto the unperturbed (and much 
denser than the gas inside the bubble) interstellar medium 
and becomes quickly radiative. 

The internal structure of the SNR is deter- 
mined by adopting th e linear velocity approximation 
|Ostriker fc McKeg 1995). in which the gas velocity profile 
for r < R sh is given by: 



i(r,t) = (l - ij u 3h (t) 



Rsh{t) 



(5) 



where a is the shock compression ratio. In the absence of 
CR acceleration, for a strong shock a = 4. Conversely, if the 
CR pressure at the shock cannot be neglected, the com- 
pression ratio increases due to the reaction of CRs onto 
the shock structure. To date, there is a general consen- 
sus on the fact that the increase of the shock compres- 
sion in SNRs is quite moderate, with a < 10, where the 
larger values are obtained for ext remely high acceler ation 
efficiencies only (e .g. Caprioli et al. 2009; Elliso n et alll201ol ; 
iBerezhko fc Volkl l20ld: IZirakashvili fc Aharonianl |201Ch . 
Here we follow IZirakashvili fc Ptuskinl 1 20121 ) and adopt 
a = 6 as a reference value. Equation [S] can be combined 
with the gas continuit y equation to determine the density 
profile inside the SNR (|Ostriker fc McKeeil 19951 ). Such pro- 
file will be used in the following to compute the hadronic 
contribution to the gamma-ray emission of SNRs. 



2.2 Cosmic ray acceleration at supernova remnant 
shocks and related gamma— ray emission 

Let f(r,p,t) be the CR particle distribution function at 
a given time t after the supernova explosion and at a 
given position r < R s h inside the SNR. Here, p is the 
particle momentum. We assume that particles are accel- 
erated at the shock with a power law spectrum so that 
/o(p,t) = f(R s h,P-,t) = A(t)p~ a , where a is treated as 
a free parameter with values in the range a > 4 (see e.g. 
IZirakasivili fc Ptuskinll2008al : ICapriolill2012l ). Values signifi- 
cantly larger than a = 4 are needed to reproduce the slope 
of the CR spectrum observed at the Earth. In the following 
we will consider two representative cases: a = 4.4 (soft spec- 
trum) and q = 4.1 (hard spectrum). The normalization A 
of the CR spectrum is computed by assuming that the CR 
pressure at the shock is equal to some fraction £cr of the 
shock ram pressure: 



p o 



- t,CR Qup U s 



(6) 



with g up representing the gas mass density upstream of the 
shock. The parameter £cr is a way to express the accelera- 
tion efficiency at the shock. For spectra steeper than a — 4 
the dominant contribution to the CR pressure comes from 
particles with momentum p fa m v c. On the other hand, the 
maximum momentum p m ax of the accelerated particles is 
determined here by the equality: 



, D(Pmax) ,r, 

Ush 



(7) 



© RAS, MNRAS 000,[THT4] 



4 P. Cristofari et al. 



where D represents the momentum dependent diffusion co- 
efficient for CRs upstream of the shock. In this way, p m ax is 
defined as the momentum for which the diffusion length of 
particles ahead of the shock Id equates some fraction £ of the 
shock radius. Particles with larger momentum are character- 
ized by larger diffusion length and are assumed to escape the 
SNR. Studies of particle acceleration and escape from shocks 
sugges t that £ ~ 0.05. ..0.1 (e.g. .Zirak ashvili fc Ptuskinl 
l2008bl . and references therein) . In the following we will adopt 
the value £ = 0.1. 

The particle diffusion coefficient at the shock depends 
on the magnetic field strength and on its turbulent struc- 
ture. If the CR acceleration efficiency is high, the intensity 
of the turbulent magnetic field at shocks is expected to be 
strongly amplified with respect to the typical value found in 
the interstellar medium, and the diffusion coefficient is sup- 
pressed accordingly. The field ampl ification m ight be due 
to a CR current driven ins tability llBelll 120041) or to reso- 
nant streaming instability (|Lagage fc CesarskvHl983T ). The 
determination of the CR diffusion coefficient in the turbu- 
lent magnetic field at shocks is still an open issue. Here we 
make the assumption that, in the presence of efficient field 
amplification, the CR diffusion coefficient is of the Bohm 
type in the amplified field. The expression for the Bohm 
diffusion coefficient is D = Rlc/3 and Rl = pc/qB is the 
particle Larmor radius. In the expressions above, c is the 
speed of light, q the elementary charge, and B the ampli- 
fied magnetic field. To estimate the value of the amplified 
field we r ely on the interpr etation of X-ray data from SNRs 
made by I York et al.l (|2005l ). Starting from the observations 
of X-ray filaments in several young SNRs, they estimated 
the magnetic field intensity just downstream of the shock 
wave. They found that a fraction £b ~ 3.5% of the shock 
ram pressure QupU 2 gh is, on average, converted into magnetic 
field. By assuming that this fraction remains constant during 
the SNR lifetime, an expression for the amplified magnetic 
field strength immediately downstream of the shock can be 
derived, and reads: 



Bdown — Bo a \l I — — 
v d 



+ 1 



(8) 



shock of the form (|Zirakashvili fc PtusMnllioil ): 



D = D B 1 + 



(10) 



where Db is the Bohm diffusion coefficient and the parame- 
ter g depends on the nature of the dominant damping mech- 
anism. We fix here g = 3, as appropriate for a Kolmogoro v 
type of non-linear damping (|Ptuskin fc Zirakashvilil [20031 ) . 
Equation[l0]is valid in both regimes u s h > Vd and u s h < Vd' 
in the former, the diffusion coefficient coincides with the 
Bohm one, while in the latter it is significantly larger than 
that due to wave damping. 

It is evident from Eq. [5] that the amplified field, being 
proportional to the shock speed, decreases with time, as long 



as u s h S> Vd- 



B d 



140 



Ush 



1000 km/s 



,V2 



(11) 



After that, field amplification becomes inefficient and the 
magnetic field downstream of the shock stays constant in 
time and is equal to aBo (i.e. only the compression of the 
field at the shock is taken into account). This fact, once 
combined with Eq. [7] implies that the maximum momen- 
tum of the protons that can be confined and accelerated at 
the shock decreases with time. As an illustrative example, 
for a SNR expanding adiabatically in an uniform medium 
Equations l8l and 1111 give, for u s h 2> Vd, B oc u s h <x t -3//5 , 
which corresponds to (Eq. [TJl p m ax oc i~ 4//5 . A quantitative 
expression for the maximum energy E max — p m axC can eas- 
ily be computed and reads: 



280 Sl[ 5 n 



1/10 C 5 TeV 



(12) 



which implies that at an age of few hundred years SNRs are 
capable to accelerate particles up to the PeV domain. In the 
opposite situation in which the magnetic field at the shock 
is not amplified and stays constant in time, the maximum 
momentum of protons exhibits a very slow decline in time. 
This happens at late times and can be described as: 

2/5 / s -1/5 

/ - IT f^l 



£si_ 

no 



where Bo ~ 5 fiG is the magnetic field in the interstellar 
medium and Vd is a velocity that defines the importance 
of wave damping i n limiting the field amplification (e.g. 
ICaprioli et al.l l2010). In Eq. [8] the term under the square 
root represents the amplification due to the CR instability 
operating upstream of the shock, while the factor a mimics 
the effect of the field compression at the shock. For shock 
velocities larger than Vd the damping of the magnetic tur- 
bulence at the shock is negligible, and thus the field can be 
effectively amplified. On the other hand, damping dominates 
for sm aller velocities. Following IZirakashvili fc Aharonian] 
l|2010h we adopt for Vd the following expression: 



f84^) 1/2 "°- 2Xl ° 8n 1/2cm / S ® 
\8Trt_BQup ) 



Vd 



which is in substantial a greement with the results of m ore so- 
phisticated studies fe.g. lPtuskin fc Zirakashvih1l2003l ). More- 
over, we adopt here a diffusion coefficient for CRs at the 



1 + 0.1 



tkyr 



30 n £5 



3/1 



TeV 



(13) 



Similar estimates can be obtained also for supernovae ex- 
ploding in a structured medium (i.e. wind plus bubble), 
though analytic expressions cannot be written in this case. 

The spatial distribution of CRs inside the remnant can 
be computed by solving the transport equation: 

df 
dt 



uVf - VDV/ - I Vui^ = 
J J 3 dp 



(14) 



where D is a momentum dependent diffusion coefficient for 
CRs. Particles with momenta smaller than p m ax are ex- 
pected to be well confined within the SNR. Thus, Equa- 
tionQ3]can be solved by dropping the diffusion term VZJV/, 
which is expected to be negligible when compared to the ad- 
vection term uVf. The solution of this differential equation 
can be found by using the method of the characteristics and 
adopting the boundary condition f(R s h,p,t) ~ fo{p,t). 

The acceleration of electrons at shocks proceeds at the 
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same rate as for protons, but their spectrum is different 
because electrons suffer synchrotron and inverse Compton 
losses, while protons are virtually loss-free. At low energies, 
where energy losses can be neglected, the acceleration of 
electrons and protons at the shock proceeds in an identical 
manner, implying that the same spectral shape is expected 
for both species. Thus, it is possible to introduce a param- 
eter K ep , generally believed to be much smaller than unity, 
that describes the ratio between the electron and proton 
spectra at low energies. 

The maximum energy of the electrons accelerated at 
a shock can be obtained by equating the acceleration rate 
at the shock to the synchrotron energy loss time. To com- 
pute t his, we follow the approach described in lVannoni et al.l 
(2009), which gives: 



Ush 



Bdown 

fOOO km/s J V 100 fiG 



-1/2 



TeV (15) 



At this energy, a cutoff appears in the electron spectrum, 
with shape oc exp [— {E/E^ nax ) 2 ] (jZirakashvili fc Aharoniar] 
l2007h . 

Electrons are accelerated very quickly, over time scales 
significantly shorter than the synchrotron energy loss time 
which, for 10 TeV electrons in a 100 piG field is of the or- 
der of a century. After being accelerated they are advected 
downstream of the shock, where they continue to lose enegy 
mainly through synchrotron radiation, with a characteristic 
time: 



1.8 x 10' 



(tbv) 



^dou 



100 /iG 



vi- 



ae) 



where E e is the electron energy. The energy loss time de- 
creases with particle energy and this implies that an energy 
Ebreak exists above which the loss time is shorter than the 
SNR age T age - Above such energy, the electron spectrum 
is shaped by radiati ve losses and steepens by one power in 
momentum (see e.g. lMorlino fc Cap rioli 2012J). In fact, since 
the SNR is expanding also adiabatic losses have to be taken 
into account, wit h a rate r a d = R s h/u s h- After a comparison 
with the work of lFinke fc Dermerl (|2012l ) we found that an 
appropriate expression for E^ reak can be found by solving 
the equation T~ g \ = r~y n + T~j. 

Two things have to be noted. First, if the magnetic field 
is not strong enough, E^ reak can easily become larger than 
Emax- In this case, no break appears and the electron spec- 
trum has the same shape as the one of protons up to Ef nax . 
Secondly, in some situations Eq. [7] can be more stringent 
than Eq. [15] (i.e. the acceleration of electrons is limited by 
escape rather that by energy losses), and in this case the 
former is used to estimate the energy of the cutoff in the 
electron spectrum. 

The last missing piece of information is the value of the 
parameter £ck, defined in Eq.JfJ which represents the parti- 
cle acceleration efficiency at the shock. In order to estimate 
this parameter, we assume that SNRs are the sources of 
galactic CRs. The estimated CR luminosity of the Galaxy is 
of the order of Lqr w 10 41 erg/s and this number is quite 
stable with respec t to the assumptio ns made to derive it 
jDogiel et al.ll2002l : IStrong et alJ^OlCf Fl By combining this 



information with the estimat ed rate of sup ernovae in the 
Galaxy vsn ~ 3/century (e.g. iLi et al.ll201ll . and references 
therein) it is possible to obtain the average fraction of the 
supernova explosion energy that needs to be converted into 
CRs in order to provide the required power Lqr . This gives: 



VCR = 



^cr v sn 

£,5N 



0.1 



/ tMW 
I ^CR 



VSN \ / £SN 



Vl0 41 ergy \0.03yr- 1 J Vl0 51 erg 

(17) 

The parameter t]cr represents the global (i.e. integrated 
over the whole SNR lifetime) CR output from a single SNR, 
while the parameter £cr that appears in Eq. [6] measures 
the instantaneous (i.e. at a specific time) acceleration effi- 
ciency at the shock. I n spired by t he result s pres ented by 
IZirakashvili fc Ptuskinl (|2012h and ICapriolil i|2012l ) we will 
assume in the following that ^cr remains constant over the 
SNR lifetime up to the end of the Sedov phase and that 
£cr ~ r\CR ~ 0.1. However, the value of the CR accelera- 
tion efficiency which is expected from theoretical studies of 
non-linear shock acceleration is generally significantly larger 
that the modest t)cr ~ 10% needed to sustain the observed 
flux of galactic CRs. One way to solve this apparent discrep- 
ancy is to assume that CRs can be accelerated with high 
efficiency (and thus modify the shock str ucture) only in a 
small fraction of the shock surface (see e.g.lVolk et al.ll2~003l ; 
iBerezhko et~ai] 120091 : IZirakashvili fc Ptuskinl l2012h . Thus, 
here the value £cr ~ 0.1 has to be considered as an av- 
erage over the whole SNR shock surface. 

The procedure described in this section allows to deter- 
mine the spectrum and the spatial distribution of CRs inside 
a SNR. By combining these results with the spatial distri- 
bution of the gas inside the SNR as obtained in Sec. 12.11 
the gamma-ray luminosity from a given SNR can be com- 
puted, by adding t he hadronic contrib ution from proton- 
proton interactions (K elner et al.| [2006f) to the leptonic one 
from inverse Compton s cattering off photons in the cosmic 
microwave background dBlumenthal fc Gould! fl970l ). While 
computing the gamma-ray emission from proton -proton in- 
teractions the results from iKelner et al.l (|2006h have been 
multiplied by a factor of ~ 1.8 to take into account the 
presence of n uclei heavie r than Hydrogen in both ambient 
gas and CRs l|Morill2009l '). 



3 ON THE NUMBER OF SUPERNOVA 

REMNANTS DETECTABLE IN THE TEV 
ENERGY DOMAIN: A MONTE CARLO 
APPROACH 

In this section we describe the simulation procedure adopted 
to predict the number of galactic SNRs with a given gamma- 
ray flux. Since this work is focused on the TeV energy do- 
main in the following we will always consider integral fluxes 
computed above a photon energy of 1 TeV. 

A Monte Carlo approach is used to simulate the time 
of explosion of all the supernovae in the Galaxy (i.e. the 
age of all the SNRs in the Galaxy). This has been done 
by assuming a supernova explosion rate constant in time 



10 40 erg for CRs in the energy range 0.1... 100 GeV. IPogiel et all 



1 IStrong et all j2010l ) give luminosities in the interval 6. 6. ..8.1 X j2002f ) give values of in the range 0.5...1 X 10 41 erg. 
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Type 


£51 


M ej , e 


M-s 


H"w ,6 


Rel. rate 


la 


1 


1.4 






0.32 


IIP 


1 


8 


1 


1 


0.44 


Ib/c 


1 


2 


1 


1 


0.22 


lib 


3 


1 


10 


1 


0.02 



Table 1. Supernova parameters adopted in the simulation: su- 



,. rg 



pernova type (column 1), explosion energy in units of 10 51 
(column 2), mass of ejecta in solar masses (column 3), the wind 
mass loss rate in M©/yr (column 4), the wind speed in units of 
10 km/s (co lumn 5), and the rel ative explosion rate (column 6). 
Values from lPtuskin et alj <201C|) . 



and equal to usn = 3/ century (|Li et al.l 1201 il l . This im- 
plies an acceleration efficiency at the shock of the order of 
rjcR ~ 10% (see Eq. ll7[) . Different values of vsn will be also 
discussed in the following. Once the age of a SNR is known, 
its location within the Ga laxy is determined by following th e 
prescription described in IFaucher-Giguere fc KaspH |2006). 
This consists in assuming that the radial distribution of 
SNRs in the Gala xy follows the distribution of pulsars , 
as determined by (|Yusifov fc Kiicukl |2004l : iLorimer] I200I ) . 
In addition to that, four spiral arms are considered, each 
arm following a logarithmic sp iral shape (see Table 2 in 
IFaucher-Giguere fc Kasp i 2006). T he wid th of the arms has 
been modeled as in iBlasi fc Amatd (|2012l ). To determine the 
height above (or below) the galactic plane of a SNR, we as- 
sume that the vertical distribution of supernovae follows the 
one of the gas. We use the vertical distribution of molecular 
and atomic Hydrog en for core-collapse and type la super- 
novae, respectively (|Shibata et al.ll201ll ). which implies that 
the distribution of type la supernovae has a height above 
the disk which is significantly larger than the one of core- 
collapse supernovae. In the absence of a detailed knowledge 
of the spatial distribution of supernovae of a given type, 
this assumption mimics the fact that core-collapse super- 
novae are expected to explode in dense star forming regions, 
while thermonuclear ones can be found also in low density 
regions. 

The dynamical evolution of each simulated SNR is then 
determined as explained in Sec. 12.11 The evolution depends 
mainly on the value of the ambient density at the location 
of the SNR and on the supernova type. To determine the 
value of the ambient density we use the three-dimensional 
distributions (galactic latitude, longitude, and radial veloc- 
ity) of atomic (HI) and molecu lar (H2) Hydrogen given by 
iNakanishi fc Sofud (|2003l . l2006h . The three-dimensional spa- 
tial distribution of the gas (i.e. the c onversion from radial ve - 
locity to distance) is computed as in lCasanova et all l|2010h . 
Four types of supernovae are considered: la, IIP, Ib/c, and 
lib, wi th relative rates 0.3 2, 0.44, 0.22, and 0.02, respec- 
tively (|Ptuskin et alj|20ich . The parameters used for each 
supernova type to compute the SNR dynamical evolution 
are listed in Table [T] The gamma-ray emission from each 
SNR is then computed as described in Sec. 12.21 

The procedure described in this Section can be used 
to simulate the number of SNRs that one can expect to 
detect with a Cherenkov telescope with a given sensitivity. 
This will be done in the next Section, where the prediction 
from our Monte Carlo will be compared with the data from 
the survey of the galactic plane performed by the H.E.S.S. 



collaboration. Before doing that, we compute here the total 
number of SNRs in the Galaxy which are expected to emit 
gamma-rays above a given flux. All the results reported in 
the following have been obtained by averaging 1000 Monte 
Carlo realizations of the Galaxy. 

We first consider a situation in which the magnetic field 
at the shock is not amplified. In this case, particles are ac- 
celerated at a shock characterized by an upstream magnetic 
field of B up — Bo « 5 fiG and a downstream magnetic field 
of Bdown = &Bo ~ 30 fiG. In Figure[T]we plot the number of 
SNRs in the Galaxy which are expected to have an integral 
gamma-ray flux above a given value. Integral fluxes above 
1 TeV are considered. The red solid line shows our prediction 
for a soft spectrum of accelerated CRs with slope a = 4.4. A 
very small electron-to-proton ratio K ep — 10~ 5 is assumed, 
and this insures that for all the SNRs the hadronic emission 
largely dominates over the leptonic one. The shaded red re- 
gion around the curve shows the fluctuations of the results 
due to the stochasticity of the process. To estimate that, his- 
tograms representing the number of realizations that corre- 
spond to a given number of detections have been produced 
and fitted with continuous functions. The shaded region rep- 
resents the interval within which 68.2% of the area below the 
fitting function is contained. The black dashed curve (as well 
as the shaded black region) has been computed, instead, by 
assuming a high electron-to-proton ratio K ep — 10~ 2 . The 
number of SNRs expected at each flux is increased by a fac- 
tor of > 1.5 with respect to the case K ep = 10 -5 , and this 
is due to the fact that the leptonic contribution is no longer 
negligible. The fact that the increase in the number of TeV- 
bright SNRs is modest but not negligible (it is indeed close 
to a factor of 2) indicates that for K ep = 10 -2 the num- 
ber of SNRs for which the hadronic emission dominates the 
gamma-ray flux is of the same order of the number of SNRs 
for which the leptonic emission dominates. 

The number of SNRs with an integral flux greater than 
1% of the Crab (« 2.3 x 10~ 13 cm~ 2 s _1 , lAharonian et all 
l2006d ), a representative flux sensitivity for deep pointed 
observations with Cherenkov telescopes of current genera- 
tiorfl is ~ 13 and ~ 21 for the red and black curve, respec- 
tively. On the other hand, the probability to detect very 
bright SNRs with fluxes of the order of ~ 10 -11 cm -2 s _1 
is small, but not completely negligible. For this value of 
the integral flux, the mean values for the expected num- 
ber of gamma-ray SNRs are ~ 0.2 and ~ 0.5 for the 
red and black curve, respectively, while the shaded regions 
(representing one standard deviation) extends up to « 1.3 
and « 1.7. This is still consistent, though in some tension 
with the fact that two SNRs have been detected at such 
flux levels: RX J1713. 7-3946, with an in tegral flux equal 
to F(> 1 TeV) » 1.6 x 10" 11 cm" 2 s _1 ijAharonian et all 



l2006bl ) and Vela Jr, with an inte gral flux equal to F(> 

1 TeV) w 1.5 x 10~ n cm" 2 s _1 ijAharonian et al.ll2007n . 
Unfortunately, a more quantitative comparison between our 
predictions and available data is, at this stage, not easy to 

2 This is true as long as point sources are considered. For ex- 
tended sources, as SNR often are, the sensitivity is worse, and 
roughly scales as the source size. So, the numbers reported here 
have to be considered only as reference (and quite optimistic!) 
values. A full discussion of this issue can be found in Section [4] 
below, where a comparison with existing data is attempted. 
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Figure 1. Number of SNRs in the Galaxy with integral gamma- 
ray flux above F(> 1 TeV). The spectral slope of accelerated 
particles at the shock is a = 4.4, the electron— to-proton ratio is 
K ev = 10 — 5 and 10 -2 for the red and black curve, respectively 
(models Ml and M2 in Table EJ). 

be performed, due to the lack of a complete catalogue of TeV 
gamma-ray bright SNRs. This point will be extensively dis- 
cussed in Section [4] 

A much m ore plausible scenar io, supporte d by both the- 
ory |Bell 2004) and observations (jVink 2012), is the one in 
which the magnetic field at the shock is substantially am- 
plified due to CR-induced instabilities that may operate in 
the shock precursor. In this case the values of the magnetic 
field and of the maximum energy of accelerated protons 
can reach values up to hundreds of microGauss and PeV 
energies or even more, respectively. These high values are 
achieved early in the evolution of the SNR and then gradu- 
ally decrease with time as the shock slows down. A plausible 
parametrization of this behavior has been described in Sec- 
tion 12.21 The expected number of SNRs with integral flux 
above F(> 1 TeV) is shown in Figure for K ep = 10" 5 
and 10 -2 (red and black curve, respectively). The dashed 
regions have the same meaning as in Figure [T] The number 
of SNRs with flux above the 1% of the Crab is ~ 32 and 
~ 39 for the red and black curve, respectively. The mean 
value for the number of very bright SNRs, with fluxes above 
10 _11 cm _2 s _1 is, for both curves, ~ 1, in closer agreement 
with the detection of the two very bright SNRs. 

The first thing to be noted is that for K ep = 10 _J (i.e. 
no leptonic contribution to the gamma-ray emission) the 
number of SNRs at a given flux is larger (by roughly a fac- 
tor of ~ 2. ..4) when the magnetic field is amplified. This 
can be seen by comparing the red lines in Figures [1] and [2] 
The reason for this is the fact that, in order to detect the 
hadronic interaction of a SNR above a photon energy of 
1 TeV, the underlying proton spectrum must extend up to 
energies significantly larger than ~ 10 TeV, because these 
are the particles that produce the photons with energy in 
excess of 1 TeV. The maximum energy of the protons accel- 
erated at the shock is larger if the field is amplified and thus, 
in this case, SNRs remain visible above 1 TeV for a longer 
time. This explains why one expects to see more gamma-ray 
SNRs if the magnetic field is amplified (even if the acceler- 
ation efficiency is the same in the two cases). 



Figure 2. Same as in Fig. [T] with the only exception that an am- 
plified field has been considered (models M3 and M4 in Table [2{. 



Another thing to be noted is that, if the field is ampli- 
fied, there is not much difference in our predictions if a low or 
a high value of the electron-to-proton ratio K ep is adopted. 
In fact, the black and red curves in Figure [2] which refer to 
K ep = 10~ 2 and 10 -5 , are virtually identical for gamma-ray 
fluxes larger than > 10 _12 cm _2 s _1 , and remain comparable 
for all the values of the gamma-ray fluxes (the difference 
between the two curves is always less than w 30%). This 
implies that the leptonic gamma-ray emission from SNRs 
never plays a crucial role. If the field is amplified, electrons 
suffer severe synchrotron energy losses, and as a consequence 
of that, their spectrum exhibit a break at an energy which 
can be computed by equating the energy loss time with the 
SNR age (see the discussion following Eq. ll6p . For large val- 
ues of the field (few hundreds microGauss) and typical SNR 
ages of thousands of years the break appears at TeV ener- 
gies. The electron spectrum below the break is identical to 
the proton one (i.e. it is a power law in momentum with 
slope 4.4) while above the break the spectrum steepens by 
one power in momentum. Such steepening suppresses the 
leptonic emission in the TeV domain, and explains why the 
parameter K ep virtually plays no role in this case. 

In computing the curves in Figure [2] we assumed that 
the magnetic field downstream of the shock is the amplified 
one, as determined by Eq. 1111 This might not be the case if 
the turbul ent magnetic fiel d is damped downstream of the 
shock (e.g. lPohl etai1l2005f ). This fact led lAtovan fc Dermerl 
i|2012l) to build a two-zone model for SNRs in which parti- 
cles are accelerated in a small region around the shock wave 
(zone 1), where the magnetic field is amplified. Particles are 
subsequently transported (through advection and diffusion) 
further inside the SNR (zone 2), where the magnetic field 
strength may be smaller. In Atoyan & Dermer's model, the 
acceleration region is much smaller than the inner region, 
and thus the electrons spend most of the time in the latter. 
Interestingly, this allowed them to decouple the region where 
electrons are accelerated (zone 1, where the magnetic field 
strength is large) from the one in which they suffer most of 
the synchrotron losses (zone 2, where the field strength is 
smaller). Assuming the existence of two zones with differ- 
ent magnetic field can significantly affect the predictions of 
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Figure 3. Same as in Fig.[TJand[2]but now two zones are consid- 
ered: the acceleration zone around the shock, where the magnetic 
field is amplified and given by Eq. 1 1 11 and an inner region with a 
weaker magnetic field equal to 30 fj,G. The spectrum of particles 
accelerated at the shock is a power law in momentum with slope 
4.4 and 4.1 for the black and red curve, respectively (models M5 
and M6 in Table fj}. 



the electron spectrum and thus of the leptonic gamma-ray 
emission from SNRs. Here we adopt the following simplified 
view: we neglect particle diffusion and we assume that elec- 
trons, after being accelerated in zone 1 are quickly advected 
into zone 2, which is characterized by a low field. In this 
case the maximum energy of accelerated electrons is com- 
puted through Eq. 1151 where Bdown is the amplified field, 
while in order to compute the energy of the break in the 
spectrum (see Eq. [16] and following discussion) we adopt a 
smaller value for the magnetic field. As an illustrative ex- 
ample, we adopt here a constant value of 30 fiG for the field 
strength in the inner region (the effect of changing this field 
will be discussed in Section [5jl. 

Results for this two-zone model are shown in Fig. [3] 
The black curve has been computed for a CR spectrum at 
injection with slope a — 4.4 and for K ep = 10 -2 . The mag- 
netic field at the shock is the amplified one (Eq.lll[) while the 
field in zone 2 is 30 fj,G. Due to the lower value of the mag- 
netic field in zone 2, synchrotron losses are less severe and 
the break in the electron spectrum moves upward in energy, 
thus enhancing the inverse Compton emission. This explains 
the larger number of gamma-ray bright SNRs expected in 
this case, when compared to the one zone model illustrated 
in Fig. [2] The number of SNRs with integral gamma-ray 
flux above 1% of the Crab flux is ~ 57 while the mean value 
for the expected number of very bright SNRs with integral 
flux above F(> 1 TeV) = Hr u cm _2 s _1 is » 1.6. 

Finally, the red curve in Fig. [3] shows the expectations 
for a hard spectrum of the accelerated CRs with a = 4.1 
(all the other parameters are left unchanged). The evident 
effect of an hard spectrum is a large increase of the num- 
ber of gamma-ray SNRs. In this case, the expected number 
of SNRs with integral flux above 1% of the Crab is un- 
reasonably large w 190, while the mean number of very 
bright SNRs with flux above 10 -11 cm _2 s -1 is « 8.1, also 
exceedingly large. This clearly disfavor a scenario in which 
SNRs accelerate hard spectra of particles. Spectra signifi- 



Model 


a 




amplified B 


# of zones 


Ml 


4.4 


10~ 5 


OFF 


1 


M2 


4.4 


io~ 2 


OFF 


1 


M3 


4.4 


10" 5 


ON 


1 


M4 


4.4 


10" 2 


ON 


1 


M5 


4.4 


10~ 2 


ON 


2 


M6 


4.1 


10" 2 


ON 


2 



Table 2. Values of the parameters adopted to compute the curves 
in Figures [TJ (model Ml and M2), [2] (M3 and M4), and E] (M5 
and M6). a is the slope of the spectrum of CRs accelerated at 
the shock, and K ep is the electron— to-proton ratio. The last two 
columns specify whether or not magnetic field amplification has 
been taken into account, and the number of zones adopted to 
compute the inverse Compton radiation from electrons (see text 
for more details). 



cantly steeper than a = 4 are needed to be consistent with 
observations, if a standard ~ 10% CR acceleration efficiency 
is assumed. This point will be further discussed in the next 
Section. 

Finally, it is instructive to estimate the total number 
of SNRs which are currently in the Sedov stage of their dy- 
namical evolution. This would provide a strict (and clearly 
over-optimistic) upper limit for the number of possible de- 
tections in gamma-rays, since CR production is believed to 
be efficient during this phase of the SNR evolution. By as- 
suming a duration of the Sedov phase equal to a few 10 4 yr 
and 3 supernova explosions per century in the Galaxy, this 
number turns out to be w 1000. Thus, for the cases con- 
sidered above, even for the most optimistic, the SNRs with 
TeV gamma-ray fluxes above the level of 1% of the Crab are 
a small fraction (~ 0.01. ..0.1) of the total number of SNR 
which are believed to accelerate CRs in the Galaxy. 

For the reader's convenience, the parameters which have 
been used to compute the curves in Figures [1] [2] and [3] are 
listed in Table fj 



4 SUPERNOVA REMNANTS AND SKY 

SURVEYS IN THE TEV ENERGY DOMAIN 

In this section we perform a comparison between the pre- 
dictions described above and the data currently available in 
the TeV gamma-ray domain. With this respect, the data ob- 
tained by the H.E.S.S. array of Cherenkov telescopes seem 
to be the most appropriate. Due to the large instrumental 
field of view (« 5°) it has been possible to devote a sig- 
nificant fraction of the total observing time of H.E.S.S. to 
a scan of the galactic plane. The re sults of this scan have 
been published in a seri es of papers l|Aharonian et al ] |2005l . 
2006a; iGast et al,ll2012h . The aim of the scan is to obtain 
a good compromise between the fraction of the sky cov- 
ered by the survey and the depth and spatial homogene- 
ity of the exposure. The original H.E.S.S. survey covered 
the rang e of \l\ < 30° in galact ic longitude and |6| < 3° in 
latitude l|Aharonian et aLll2005T ), and it has been gradually 
extended thereafter, especially in longitude. To date, an ex- 
tension in the ran ge of I = 250 to 65 degrees was reported 
|Gast et al.ll2012T ). However, from Fig. 2 in that paper, it 
can be noticed that the exposure, and thus the sensitivity 
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within the survey region is non-uniform. For this reason, 
in the following we restrict our attention to the region of 
galactic longitude |2| < 40 degrees only, within which the 
sensitivity for point sources is quite homogeneous and al- 
ways at the level of at least ~ 1.5% of the Crab level (i.e. 
F(> 1 TeV) « 3.4 x 10 _13 cm _2 s _:L ). 

The number of TeV gamma-ray sources in the H.E.S.S. 
Source Catalog within the region we selected < 40°, 
\b\ < 3°) and with a flux above 1.5% of the Crab is 35. No- 
tably, three of them, RX J1713.7-3946, HESS J1731-347, 
and CTB 37B, are associated with SNR shells. The physical 
properties of these three sources are listed in Tabel[3] In ad- 
dition to that, other three sources, CTB 37A, HESS J1745- 
303, and W28, are or might be associated with SNR shells 
in interaction with massive molecular clouds. However, the 
gamma-ray emission from these interacting systems might 
have a different origin than the one we investigate here. For 
example, the gamma-ray emission from the old SNR W28 
(the estimated age is a few times 10 4 yr) has been interpreted 
as the results of the interactions of CRs that escaped the 
SNR and penetrated into t he molecular cloud l|Gabici et al.l 
l20ld ; iNava fc Gabidl20l3 ). Also the SNR coincident with 
the gamma-ray source HESS J1745-303 is believed to be 
quite old (more than ~ 10 4 vr. lAharonian et al.ll2008bT ). and 
thus also in this case an interpretat ion of the gamma- ray 
emission in terms of escaping CRs (|Gabici et alj|2009h or 
of re-accelerati on of pre-existing CRs (the so called cloud 
crushed model, lUchivama et alj |2010| ) might be preferred. 
The situation is different for the S NR CTB 37A, which is 
likely to be young (« 1...3 x 10 3 vr. lAharonian et al.ll2008d 1 
and thus in this case the gamma-ray emission might be re- 
lated to the ongoing acceleration of CRs at the SNR shock, 
though the presence of the cloud might significantly affect 
the general picture of diffusive shock acceleration described 
in Sec. 12.21 However, the gamma-ray emission from CTB 
37A can also be attributed to a pulsar wind nebula present 
in the region. For these reasons we do not consider these 
three objects in the following. Finally, 17 out of the 35 TeV 
sources detected by H.E.S.S. in the region in exam still re- 
main unidentifiecy. Thus, we can conclude that the number 
of SNRs detected in TeV gamma rays in the considered re- 
gion spans from a pessimistic tally of 3 (if only the three 
isolated shells listed in Table [3] are considered) to an overop- 
timistic one of ~ 20 (in the unlikely event that most or all 
of the unidentified H.E.S.S. sources are indeed SNRs). 

To compare these numbers with our predictions we run 
1000 Monte Carlo realizations of the supernova explosions 
in the Galaxy and compute the number of sources expected 
to be detected by H.E.S.S. within the region in exam. A 
sensitivity at the level of 1.5% of the Crab flux has been 
adopted for point like sources, while for extended ones the 
sensitivity has been degraded by a factor of "ds/fipsF, where 
•& s is the source apparent size and $psf ~ 0.1° is the angu- 
lar resolution of H.E.S.^f). The results of this computation 



3 www.mpi-hd.mpg. de/hf m/HESS/pages/home/sources/ 

4 The classification of TeV sources in shells, shells interacting 
with molecular clouds, and unidentified sources has been estab- 
lished by cross— correlating the H.E.S.S. source catalogue with the 
TeVCat online catalogue, maintained by S. Wakely and D. Horan 
http : //tevcat .uchicago . edu/ 

5 A discussion of the procedure to determine the extension of a 
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Figure 4. The top panel shows the expected number of SNRs 
detectable by H.E.S.S. in the region of coordinates \l\ < 40° and 
|b| < 3° as a function of the spectral slope a of accelerated par- 
ticles. Black and red lines correspond to K ep = 10~ 5 and 10 -2 , 
respectively. The other panels (top to bottom) show, for the case 
K ep = 10 — 2 the distance, age, and angular size of the detected 
SNRs. Median values of these quantities are shown for all SNRs 
and for spatially resolved ones with a solid and dashed thick line, 
respectively. The thin dashed lines represent the maximum value 
for these quantities (averaged over the Monte Carlo realizations). 
In the bottom panel the fraction of point-like sources is also 
shown as a black dotted line. 



are shown in the top panel of Fig. 3J where the mean num- 
ber of expected detections is plotted as a function of the 
spectral slope a of accelerated CRs. The black and red lines 
refer to values of the electron-to-proton ratio of K ep = 10~ 5 
and 10~ 2 , respectively. A two-zone model has been adopted 
to describe the SNR, with a small region around the shock 
where the magnetic field is strongly amplified (see Eq. [8]), 
and an inner region inside the SNR where the magnetic field 
is significantly lower due to damping. We assume a field in- 
tensity of 30/iG in the inner region, but we discuss in the 
next section which effect a change in the value of this pa- 
rameter has. 

It can be seen from Fig. [4] that for steep particle spec- 
tra (a — 4.4) the expected number of detections is « 5 and 
7 for K ep = 10~ 5 and 10~ 2 , respectively. These numbers 
progressively increase if harder and harder spectra are con- 
sidered, and for a — 4.1 we found, for the two values of K ep , 



TeV source clearly goes beyond the scope of this paper. However, 
it is important to remind that the classification of a sources as 
extended may depend on several factors, including the available 
photon statistics. The value 0.1° adopted here must be considered 
as an indicative figure only. 
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Name 


F(> 1 TeV) 


d 


age 


radius 


Ref. 




[10" 12 cm- 2 s- 1 ] 


[kpc] 


[kyr] 


[°] 




RX J1713.7-3946 


15.5 


1 


1.6 


0.65 


1,2,3 


HESS J1731-347 


6.9 


2.4.. .4 


27 


0.25 


4,5 


CTB 37B 


0,1 


13.2 


0.3... 3 


0.03 


6.7 



Table 3. Gamma-ray fluxes, distances, ages and apparent size s of the three SNR shell s d etected by H.E.S.S. in th e region / < 40° , 
|6| < 3.5° at a flux level ab ov e 1.5% of the Cra b. References: 11 lAharonian et al.ll2006bl; 2) iMoriguchi et~al"1l2005l; 31 IWane et aTlll997l ; 
4) lAbramowski et al.|[201ll ; 5) iTian et al.ll200Sl ; 6) lAharonian et alj|2008al ; 71 iNakamura et al.ll2009l 



~ 17 and 22 detections. These numbers are comparable, or 
even larger than the most optimistic estimate for the actual 
number of SNRs detected by H.E.S.S. in the region. This 
fact suggests that steep spectra are consistent with observa- 
tions, while hard spectra, with a close to 4 (the standard 
prediction of test-particle first order Fermi acceleration) 
would imply an unreasonably high number of detections. 
This finding is in agreement with both gamma -ray observa- 
tions of some individual SNRs like Tycho (e.g. lAcciari et al.l 
l201ll ; iMorlino fc Capriolill2012|l and recent theoretical devel- 
opments (e.g. IZirakasivili fc Ptuskmll2008al ; iBell et al.|[201ll ; 
Caprioli 2012) that seem to indicate that diffusive shock ac- 
celeration might produce spectra significantly steeper than 
a = 4 if diffusive shock acceleration is treated in a fully 
self-consistent manner. 

All the other panels of Fig. [4] refer to the case K ep = 
10 -2 , and show the median and maximum values of the fol- 
lowing quantities (top to bottom): distance, age, and appar- 
ent angular size of the SNRs that should have been detected 
by H.E.S.S. in the region we are considering. Thick red lines 
show the median values for all those SNRs (dashed lines) 
and for spatially resolved ones (solid lines), i.e., with an an- 
gular size larger than ~ 0.1°. The thin dashed lines show 
the maximum value for these quantities, averaged over the 
number of Monte Carlo realizations. 

The median distance of detected SNRs lays, for all the 
values of the spectral slope a, in the range 5 kpc < d < 
10 kpc, which means that in the majority of the cases the 
detected SNRs are closer than the galactic centre. The me- 
dian distance is slightly smaller (d ~ 5 kpc) if only resolved 
sources are considered, as expected given the worse instru- 
ment sensitivity in detecting extended sources, and given 
that it is easier to resolve nearby sources. On the other hand, 
the maximum distance up to which SNRs are detected - a 
sort of horizon for the detection of SNRs - is ~ 15 kpc for 
hard particle spectra (a ~ 4) and decreases gradually for 
steeper and steeper spectra reaching a value of « 10 kpc for 
q = 4.4. 

The predicted median age of the SNRs detectable by 
H.E.S.S. is quite insensitive to the value of a, and is of the 
order of < 5 kyr. This slightly increase to > 5 kyr if only 
resolved SNRs are considered. This is expected, given that 
older SNRs are obviously larger than younger ones. Also in 
this case, the maximum age of the detected SNRs is pre- 
dicted to decrease from m 20 kyr to « 12 kyr when the 
spectral slope of accelerated CRs goes from a = 4.1 to 4.4. 

Finally, the predicted fraction of point-like SNRs is in 
the range « 0.4. ..0.6, as indicated by the black dotted line in 
the bottom panel of Fig. [4] Amongst extended sources, the 
expected median angular size is ~ 0.2°, while the largest de- 



tectable sources have a size of « 1...1.2". All these quantities 
are quite insensitive to the value chosen for a. 

Though a rigorous comparison between our predictions 
and available data is not easy, it is evident that a qualitative 
agreement between data and predictions exists. Our expec- 
tations for the selected region of the H.E.S.S. scan seem to 
reproduce correctly the actual number of detections and the 
typical distances, ages, and apparent sizes of the gamma- 
ray bright SNRs (see e.g. Fig. [4] and Table [3}. Moreover, as 
already discussed in Sec. O also the number of very bright 
(flux at the level of « 10 -11 cm~ 2 s _1 ) SNRs detectable in 
the whole Galaxy seems to be well reproduced. All these 
facts are encouraging and provide additional support to the 
consistency of the SNR paradigm for the origin of CRs. A 
summary of the main finding for the different scenarios con- 
sidered in this paper can be found in Table [4] 

We conclude the Section with two more predictions 
of our calculations. The first one concerns the fraction of 
gamma-ray bright SNRs whose emission is dominated by 
hadronic processes. This fraction, as can be seen from Ta- 
ble 21 depends quite strongly on the adopted parameters 
(especially on the magnetic field strength). It can range 
from m 60% to 100%. Determining the hadronic or leptonic 
origin of the gamma-ray emission from a given SNR is a 
very difficult task. Consider, for example, the three SNR 
listed in Table [3] While multi-wavelength observations of 
RX J1713. 7-3946 seem to point towards a leptonic origin of 
the gamma-ray emission (but see lFukui et al . 2012 for an al- 
ternative explanation), for the other two SNRs the situation 
is still ambiguous. Thus, if the commonly accepted interpre- 
tation of the observations of RX J1713. 7-3946 is correct, 
at least for some SNRs the detected gamma-ray emission 
must be leptonic, and this would disfavor models Ml and 
M3, for which the electron-to-proton ratio is very small 
(K ep — 10 ), and also model M4 for which a very large 
magnetic field has been assumed. The second prediction is 
the fact that a very large fraction (« 80% for model M5, 
~ 65% for model M2) of the SNRs which are expected to 
be detected in TeV gamma-rays are of type la. The dif- 
ficulty of detecting core-collapse supernovae is connected 
to the fact that the SNR shock propagates in the tenuous 
medium of the wind-blown bubble, which strongly reduces 
the gamma-ray emission due to proton-proton interactions 
(the core-collapse SNRs which are expected to be detected 
in gamma-rays are characterized by a dominant leptonic 
emission). Determining the type of the progenitor supernova 
is a very difficult task. Amongst TeV-bright SNRs, only a 
few have been firmly identified as thermonuclear or core- 
collapse supernovae. From the detection of the light echoes 
of the supernova explosions, Tycho has been firmly identified 
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Model: 


Ml 


M2 


M3 


M4 


M5 


M6 


Mean (median) number of detections: 


0.9 (2) 


1.8 (3) 


5.3 (6) 


5.9 (6) 


6.6 (7) 


22 (23) 


Median distance [kpc]: 


2.6 


2.7 


5.0 


5.3 


5.0 


8.7 


Median age [kyr]: 


1.8 


1.0 


4.2 


3.0 


2.8 


4.2 


Median apparent size [°]*: 


0.25 


0.28 


0.22 


0.26 


0.22 


0.20 


Fraction of point sources: 


0.34 


0.41 


0.40 


0.51 


0.40 


0.55 


Fraction of hadronic sources: 


1 


0.59 


1 


0.98 


0.87 


0.71 



* Extended sources only (i.e. size larger than 0.1°). 



Table 4. Expected number of detections in the region considered in Fig.|4]and SNR properties for the different models listed in Tablc[2] 



as a type l a supernova |Krause et al.ll200ct ). while Cas A as 
a type lib (|Krause et alj|2008h . Also SN 1006 is confidently 
identified as a type la supernova d ue to its location quite 
distant from the galactic disk (e.g. IStephensonll201(il ). For 
the other TeV SNRs the situation is less clear, t hough some 
hints have been provided in some cases (see e.g. iTian et aD 
120101 ). and thus more efforts are needed in order to increase 
the number of firm identifications of the supernova type. 



5 DISCUSSION OF THE RESULTS 

In this Section we investigate how the results obtained in 
this paper change when different assumptions are made on 
the values of some key physical parameters. 

One of the most important quantities involved in our 
calculations is the global CR luminosity of the Galaxy, 
Lcr > which represents the total energy output from all the 
sources of CRs in the Galaxy. If SNRs are the main sources 
of CRs, then from the supernova rate in the Galaxy vsn 
it is possible to constrain the typical amount of energy that 
each SNR must convert into CRs (see Eq. ll7[) . The CR lumi- 
nosity of the Galaxy is determined by modeling the escape 
of CRs from the Galaxy, and different approaches lead to 
very similar values of this quantity, which is of the order 
of « 10 41 erg/s (e.g. iDogiel et al.ll2002l ; Istrong et all 

|2010| ). However, an uncertainty up to a factor of ~ 2 might 
still be accepted. For this reason, we repeat our calculations 
for three values of , namely, 5 x 10 40 erg/s, 10 41 erg/s, 
and 2x 10 41 erg/s, and, as done in Sec.[4l we compute the ex- 
pected mean number of SNRs detectable within the H.E.S.S. 
survey of the galactic disk in the galactic longitude and lat- 
itude ranges \l\ < 40° and |6| < 3°, and with integral flux 
above 1.5% of the Crab. For the three values of L^ii dis- 
cussed above the mean number of detections scales roughly 
linearly and reads 3.1, 6.6, and 11, respectively (for model 
M5 in Table [2} . The approximate linearity between the av- 
erage number of detections and the CR power in the Galaxy 
can be easily understood as follows: if we keep all the other 
parameters unchanged, the effect of varying the global CR 
luminosity is reflected into a different acceleration efficiency 
Vcr that the SNRs must achieve in order to sustain the CR 
intensity in the Galaxy. If Lqr is increased by a factor of 
/, then also the acceleration efficiency t\cr is augmented by 
the same factor, as well as the expected gamma-ray lumi- 
nosity from each SNR. This in turn implies that SNRs would 
be visible by the same telescope up to distances d a factor 
of ~ /V 2 larger, or, if as a first approximation we assume 
SNR to be homogeneously distributed in a flat disk, within 



a volume a factor of oc d 2 oc / larger, which explains the 
linearity. 

Note that an identical linear scaling has to be expected 
also if we relax the assumption of equality t)cr ~ £cr be- 
tween the two CR acceleration efficiencies defined in Sec. 12.21 
and we substitute it with the expression tjcr = f (,cr, 
where / accounts for possible deviations from the equality. 
However, as already said above, theoretical investigations 
indicates that / should be q uite close to 1 l|Caprioli 1 20121 ; 
IZirakashvili fc Ptuskinir2012h . and thus the predictions re- 
ported here can be regarded as fiducial estimates, easy to 
be rescaled for possibly different values of /. 

Another crucial parameter is the supernova explosion 
rate in the Galaxy, vsn- We adopt throughout the paper 
a value of usn = 3 /centu ry, in agreement with recent esti- 
mates (e.g. iLi et al.ll20l"ll ). However, also in this ca se a sys- 
tema tic uncertainty of a factor of ~ 2 is expected (|Li et al.l 
2011). If we repeat the estimate of the mean number of de- 
tections (model M5, |/| < 40°) for explosion rates in the 
range vsn = 1...3/century we do not obtain any significant 
difference in the predicted number of detections. This can 
be understood by noting that a change in the supernova rate 
also affects the CR acceleration efficiency per SNR. If the 
rate of supernova explosions vsn is multiplied by a factor of 
/, the acceleration efficiency per SNR (and thus its gamma- 
ray emission) has to be multiplied by the inverse factor 
in order to keep the CR luminosity in the Galaxy constant. 
In other words, there will be / more SNRs that contribute 
to the CR intensity in the Galaxy, but each SNR will be 
a factor of / less powerful in gamma-rays, and thus visi- 
ble within a volume a factor of / smaller. And this explains 
why our predictions are quite insensitive to the choice of 
the parameter vsn- This can be restated in another way: if 
we reduce vsn, the fraction of SNRs that can be detected 
by a given instrument is larger, even if the total number of 
detections is unchanged. 

Also the value of the magnetic field plays a crucial 
role and it is thus mandatory to investigate how our pre- 
dictions change if different values of the field are assumed. 
A smaller magnetic field strength reduces the synchrotron 
losses of electrons that can thus radiate more gamma-ray 
photons through inverse Compton scattering. We consid- 
ered the two-zone model (M5) and varied the intensity of 
the field in zone 2. The mean number of expected detections 
in this case goes from m 8.6 to w 6 if the field is assumed to 
vary from 5 to 40 /iG. Moreover, the fraction of SNRs whose 
TeV emission is dominated by the hadronic component goes 
from 63% to 97%. 

We investigate now which effect has on our predictions 
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a different assumption concerning the maximum energy up 
to which particles can be accelerated at SNR shocks. The 
estimates for the maximum energy given in Eq. [T2] and \T3\ 
are very plausible guesses, but it is true that we are far from 
having a solid knowledge of the details of the amplification 
mechanism of the magnetic field that determines the maxi- 
mum particle energy at a shock. A way to change the value 
of the maximum energy of accelerated particles is to change 
the size of the CR shock precursor, i.e. the value of the pa- 
rameter £ in Eq. [7] If we go from £ = 0.1 to £ = 0.05 we 
reduce by a factor of 2 both the size of the precursor and 
the value of the maximum energy. The number of expected 
detections is, in the two cases, sa 6.6 and w 5, respectively. 
Thus, we can conclude that our predictions are not much af- 
fected, unless the assumed values for the maximum particle 
energy are varied significantly (more than a factor of 2). 

Finally, we checked for the stability of our predictions 
against variations of the spatial distribution of SNRs and 
gas in the Galaxy. We repeated the procedure for a spatial 
distribution of SNRs with and without taking into account 
the presence of spiral arms, and we used the r adial distribu- 
tion of SNRs given by ICase fc Bhattacharval l) 19981 ) instead 
of the one by lLorimerl ( 20041 ) and did not found any signifi- 
cance variation in our predictions. This is connected to the 
fact that SNRs can be detected up to quite large distances, 
of the order of w 10. ..15 kpc (see Fig. [4}, and thus the effects 
of different spatial distributions are not that important. We 
also used cylindrical symmet ric templates for th e gas dis- 
tribution in the Galaxy (as in lShibata et al.ll201l[ ), without 
finding an appreciable effect. However, it has to be noted 
that the surveys of CO and HI that are used to infer the 
gas distribution in the Galaxy are characterized by a spa- 
tial resolution along the line of sight of « 50. ..100 pc. This 
might create problems, for examples, in identifying dense 
molecular clouds which have a typical size well below 100 
pc. Though this is expected to have an effect on our esti- 
mates, we know from observations that SNRs in interaction 
with molecular clouds are generally quite old, with ages o f 
the order of ~ 10 4 yr or more (e.g. lUchivama et al.l l2010). 
while the majority of the SNRs for which we predict a de- 
tectable TeV emission have an age well below 5 kyr (see 
Fig- [J). It is not clear whether such old SNRs are capable of 
accelerating particles up to energies well above « 10 TeV, 
as needed in order to produce w TeV photons. According 
to our current knowledge of the shock acceleration mecha- 
nism, which we briefly reviewed in Sec. 12.21 it seems that 
old SNRs can, at best, marginally reach these energies. This 
suggests that, with this respect, our predictions can still be 
considered solid and reliable. 



6 CONCLUSIONS 

In this paper we performed a comparison between the expec- 
tations of the SNR paradigm for the origin of galactic CRs 
and the available data in the TeV energy domain. Instead 
of proceeding on a case-by-case study of individual SNRs, 
we aimed at studying TeV-bright SNRs as a population. To 
our knowledge, this is the first time that such an approach 
is performed. 

We started by assuming that SNRs are the main sources 
of CRs, and this allowed us to estimate the typical CR ac- 



celeration efficiency per SNR. We then used a Monte Carlo 
approach to simulate the position and time of explosions of 
the SNRs in the whole Galaxy and we computed then the ex- 
pected number of SNRs that should be detected in the TeV 
domain by the present generation of Cherenkov telescopes. 
To compare our predictions with data, we selected a region of 
the galactic disk with galactic longitude \l\ < 40°, for which 
H.E.S.S. performed a scan with a roughly spatially homo- 
geneous exposure, corresponding to a sensitivity of « 1.5% 
of the Crab. Predictions seem to agree well with available 
data, thus providing an additional consistency check of the 
idea that SNRs are the sources of CRs. 

Our main findings can be summarized as follows: first 
of all, we obtained evidence for the fact that particle spec- 
tra significantly steeper than a = 4 have to be accelerated 
at SNRs, if they indeed are the sources of CRs. The reason 
for that is the fact that hard spectra (a ~ 4) would result 
in a very strong TeV luminosity and this would be incon- 
sistent with the scarce number of SNRs currently detected 
at TeV energies. Secondly, the expected fraction of gamma- 
ray bright SNRs whose emission is dominated by neutral 
pion decay strongly depends on the assumptions made on 
the strength of the magnetic field. For the range of param- 
eters investigated in Sec. [4] this fraction spans from « 60% 
to 100%, and the largest values are obtained either for a 
very low electron-to-proton ratio or for a very large mag- 
netic field strength (one of the two conditions suffices to 
increase the fraction up to « 100%). The fact that there is 
at least one SNR (namely RX J1713. 7-3946) whose gamma- 
ray emission is commonly ascribed to inverse Compton scat- 
tering might suggest that an high electron-to-proton ratio 



(of the order of K e 



10" 



characterizes the acceleration of 



particles at SNRs and that regions where the field strength 
is not too large must exist in SNRs. Finally, according to 
our predictions we should expect to detect the same num- 
ber of extended and point-like (where with point-like we 
intend sources with a size smaller than 0.1 degree) TeV- 
bright SNRs, and supernovae of type la should account for 
a large fraction of the detections (« 60. ..80%). 

Before concluding we comment on a possible exten- 
sion of the procedure described in this paper to the GeV 
energy domain, currently probed by the Fermi and Agile 
satellites. Remarkably, a constantly increasing number of 
SNRs is being detected in the GeV energy band. Several 
of these SNRs are quite old, often radiative systems that 
show clear signatures of interacti ons with massive molecu- 
lar clouds l|Uchivama et alj|2010h . For these systems, not 
considered here, the scenario of particle acceleration and 
gamma-ray production is most likely very different from 
the one considered in this paper (see e.g [Gabici et al.ll2009l ; 
iMalkov et al]|201ll ; [Uchivama et alj|20ich . Despite this fact, 
we can anyway use the procedure developed in this paper, 
and keep in mind that the estimates obtained in this way in 
the GeV domain would be very approximate, and most likely 
lower limits only (because it won't take into account the old 
interacting systems). By considering a sensitivity for Fermi 
(integrated above 1 GeV) of few 10 _9 cm _2 s _1 in the inner 
Galaxy (where most of the detections are likely to happen) 
we obtain a number of expected detections of the order of 
several tens. 

In a forthcoming publication (Cristofari et al., in prepa- 
ration) the procedure developed in this paper will be used 



© RAS, MNRAS 000,[THl4] 



Cosmic rays and gamma rays from supernova remnants 13 



to estimate the impact that the next generation Cherenkov 
Telescope Array will have on the studies of acceleration of 
CRs at SNR shocks. 
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